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Modelling magnetic flux emergence in the solar convection zone 
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ABSTRACT 



Context, Bipolar magnetic regions are formed when loops of magnetic flux emerge at the solar photosphere. Magnetic buoyancy 
plays a crucial role in this flux emergence process, particularly at larger scales. However it is not yet clear to what extent the local 
convective motions influence the evolution of rising loops of magnetic flux. 

Aims. Our aim is to investigate the flux emergence process in a simulation of granular convection. In particular we aim to determine 
the circumstances under which magnetic buoyancy enhances the flux emergence rate (which is otherwise driven solely by the convec- 
tive upflows). 

Methods. We use three-dimensional numerical simulations, solving the equations of compressible magnetohydrodynamics in a 
horizontally-periodic Cartesian domain. A horizontal magnetic flux tube is inserted into fully developed hydrodynamic convection. 
We systematically vary the initial field strength, the tube thickness, the initial entropy distribution along the tube axis and the magnetic 
Reynolds number. 

Results. Focusing upon the low magnetic Prandtl number regime (Pm < 1) at moderate magnetic Reynolds number, we find that the 
flux tube is always susceptible to convective disruption to some extent. However, stronger flux tubes tend to maintain their structure 
more effectively than weaker ones. Magnetic buoyancy does enhance the flux emergence rates in the strongest initial field cases, and 
this enhancement becomes more pronounced when we increase the width of the flux tube. This is also the case at higher magnetic 
Reynolds numbers, although the flux emergence rates are generally lower in these less dissipative simulations because the convective 
disruption of the flux tube is much more effective in these cases. These simulations seem to be relatively insensitive to the precise 
choice of initial conditions: for a given flow, the evolution of the flux tube is determined primarily by the initial magnetic field 
distribution and the magnetic Reynolds number. 

Key words. Convection - Magnetohydrodynamics (MHD) - Sun: granulation - Sun: interior 
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1. Introduction 

It is generally believed that most of the large-scale toroidal mag- 
netic field in the solar interior is stored in the stably-stratified 
layer just below the base of the convection zone, where the 
toroidal field is a mplified by the effec ts of differential rotation 
(see, for example, Ossendriiver 2003). Through the action of 
magnetic buoyancy (|Parkerlll955l) . segments of this toroidal flux 
rise into the convection zone, eventual ly emergin g to form bipo- 
lar active regions at the solar surface (Zwaan 1985). Several nu- 
merical simulations of this flux emergence process have focused 
upon the i nteractions between convection and buoyant magnetic 
flux tubes dDorch et al.l200TtlFan et al.l2003HAbbett et al.l2004t 
ISteinetalJ 20111) . whilst others have focused upon the subse- 
quent small-scale and large-scale evolution of the e merging field, 
as it r i ses into the photosphere/chromosphere (e.g .lCheung et al. 
I2007t Qbrtosa -Andreu & Moreno-Insertis 2009) and into the 
lower corona (e.g. iMagard |200U iFanl I200U lArchontis et al.l 
12004 iMartfnez-Svkora et alj "f2008). There have also been a 
number of recent reviews of flu x emergence dynamics (e.g. 
iHood et al.11201 lHArcho"ntisll2012l) . 



Since it is not yet possible to model the whole flux emer- 
gence process (from the base of the solar convection zone into 
the solar atmosphere) in a self-consistent way, it is necessary to 
consider local models with idealised initial conditions. Usually, 
the initial magnetic field distribution corresponds to a thin hor- 



izontal flux tube. This flux tube is often twisted, although the 
components of the magnetic field that are perpendicular to the 
tube axis are usually much weaker than the axial component. 
As a result, the dominant component of the Lorentz force cor- 
responds to the radial gradient of the magnetic pressure. In cal- 
culations of this type, the local gas pressure is usually modi- 
fied so that this initial flux tube is in approximate (total) pres- 
sure balance with its non-magnetic surroundings. The tube is 
then made buoyant by specifying a particular entropy distribu- 
tion (or, equivalently, a density deficit) along its axis. In some 
models, the imposed specific entropy distribution varies along 
the axis of the tube, with the result that only one part of the tube 
is buoyant. This differential buoyancy leads to the development 
of Q-like loops of magnetic flux. When convective motions are 
present, it is not necessary to impose a variable entropy distribu- 
tion along the tube in order to obtain Q-like loops. Convective 
upflows (and downflows) naturally lead to the formation of loop- 
like structures, provided that the magnetic Reynolds number is 
large enough to ensure that the magnetic field lines are (at least 
partially) advected with the flow. 

In the absence of convective motions, there is nothing to dis- 
rupt the buoyant rise of a magnetic flux tube. The situation is 
clearly more complicated wh en the tube rises through a convec- 
tive layer. In a previous studv. lFan et al] (120031) used the anelastic 
approximation to model the interaction between uniformly buoy- 
ant magnetic flux tubes and convection in the deep layers of the 
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solar convection zone. Varying the initial field strength in their 
model, they found that weak fields tended to be disrupted by 
the convective motions. In fact the peak field strength within the 
tube had to be significantly greater than the equipartition value 
(at which the field is in energy balance with the surrounding mo- 
tions) before magnetic buoyancy could play a dominant role in 
the evolution of t he flu x tube. Similar conclusions were reached 
by lAbbett et al.l d2004f) in their related study. Although we do 
not go into more details here, it is worth noting one further as- 
pect of these anelastic models. Defining the magnetic Prandtl 
number, Pm, to be the ratio of th e magnet i c Rey nolds number 
to the (fluid) Reynolds number, iFan et all d2003l) focused ex- 
clusively upon the Pm > 1 regime, at re l atively low Reynolds 
numbers. The simulations o f lAbbett et al.1 (2004) were at higher 
Reynolds number, but most were still in this Pm > 1 regime. 
In actual fact, we would expect Pm <g 1 t hroughout the solar 
convection zone (see, e.g. lOssendriiverll2003l) . However, the low 
magnetic Prandtl number regime is numerically extremely chal- 
lenging, which presumably explains why these studies did not 
focus upon this region of parameter space. 

Models of flux emergence through fully compressible con- 
vection have also been carried out. It is possible to simulate this 
proce ss on scales that are comparable to t hat of an active re- 
gion (I Cheung et al.H2010t Istein et al.ll20llt) . However, it is not 
feasible to carry out parametric surveys in calculations of this 
size, so most previous studies have focused upon flux emer- 
gence through convective layers with a relatively small number 
of granular cells. T he key finding of the compressible model of 
iDorch et"aD (|2001) is that convective disruption removes mag- 
netic flux from the tube as it rises, thus reducing the quantity 
of flux that emerges at the surface layers. In some sense, there- 
fore, convective motions reduce the efficiency of the flux emer- 
gence process. However, it i s not clear to what extent this ef- 
fect is parameter-dependent. ICheung et al.l d2007l) considered a 
model of compressible convection that included the effects of 
radiative transfer. At some instant in time, they introduced a 
long, uniformly-buoyant magnetic flux tube into the lower part 
of the domain, modifying the local gas pressure and the local 
velocity field so that the tube was (at least instantaneously) in 
equilibrium with its non-magnetic surroundings. Under the ac- 
tion of the vertical convective motions, they found that weak 
flux tubes tend to develop a sea-serpent-like form. However, as 
in the anelastic studies, magnetic buoyancy dominates the evo- 
lution of the tube only when the peak magnetic field strength 
exceeds some threshold value. At the surface these emerging 
fields can be strong enough to modify the granula tion pattern. In 
their larger scale simulations. IStein et al.l ( 201 lb adopted a dif- 
ferent approach, introducing untwisted, horizontal magnetic flux 
through the lower boundary of their convective domain, ensur- 
ing that it had the sa me entropy as the convective upflows. Like 
Cheung et al. (2007), they found that weak fields were suscep- 
tible to convective disruption. Fields that were strong enough 
to rise to the surface through the action of magnetic buoyancy 
tended to modify the surrounding flow, producing unrealistically 
hot, large granules. 

Motivated by these previous studies, the aim of this paper is 
to investigate the parametric dependence of simulations of flux 
emergence across a convectively-unstable compressible layer in 
the "low" magnetic Prandtl number regime. Most of the com- 
pressible calculations that are described above utilise artificial 
viscosities to stabilise the numerical scheme, so it is difficult 
to define a magnetic Prandtl number in these cases although, in 
practice, it is likely that Pm a; 1 in these calculations. Therefore, 
this will be the first time that the Pm < 1 regime has been studied 



systematically in this context. Given that any initial conditions in 
local models of this type are likely to be highly idealised, it is im- 
portant to determine the extent to which the evolution of the sys- 
tem is sensitive to the precise choices that are made. We there- 
fore consider several different initial configurations, varying the 
peak field strength, the width (and twist) of the tube, in addition 
to varying the initial entropy distribution within the flux tube. 
By carrying out this systematic survey we hope to identify the 
key features of the initial configuration that determine whether 
or not magnetic buoyancy contributes towards the evolution of 
the tube. The paper is structured as follows. In the next section, 
we describe the governing equations, the model parameters and 
the choice of initial conditions under consideration. In Section 
3, we present our numerical results. In the last section, we sum- 
marise our findings and discuss some of their implications for 
flux emergence in the solar convection zone. 

2. Model setup 

2.1. Governing equations 

We consider the dynamics of a plane layer of compressible, 
electrically-conducting fluid that is heated from below. The fluid 
is characterised by various parameters (all of which we assume 
to be constant in this idealised model): the thermal conductivity 
K, the shear viscosity yU, the magnetic diffusivity rj, the mag- 
netic permeability fio, and the specific heat capacities, cy and 
cp. The gas constant, R t , is defined by /?* = Cp - cy. Choosing 
a Cartesian coordinate system in which the z-axis points verti- 
cally downwards (parallel to the constant gravitational accelera- 
tion gz), this fluid occupies the region < x < 8c/, < y < 8c/ 
and < z < d, where d is the depth of the domain. Periodic 
boundary conditions are imposed in the horizontal directions. 
The upper and lower boundaries are held at fixed temperature, 
so that T = Tq at z = and T - Tq + AT at z = d. These sur- 
faces are also assumed to be impermeable and stress-free. The 
magnetic field satisfies B x = B y = at z = and z - d, which 
corresponds to a vertical magnetic field boundary condition. 

The governing equations are those of non-ideal, compress- 
ible magnetohydrodynamics, wit h an appropriate e quation of 
state for a perfect gas (see, e.g. IBushbv et alJ l2008). We non- 
dimensionalise this system (see IBushbv et al.l I2008L for more 
details), scaling all lengths by d, whilst all time-scales are 
expressed in terms of an (isothermal) acoustic travel time, 
d/(R t To)^ 2 . Hence, the fluid velocity, u, is scaled in terms of 
the unperturbed isothermal sound speed at the upper surface, 
(R t T[)) l l 2 . The temperature, T, and density, p, are scaled by Tq 
and po respectively (where p is the density at the upper surface 
in the absence of convection). We scale the magnetic field, B, by 
(piupoR t To) l l 2 . This choice of scaling for the magnetic field im- 
plies that the Alfven speed at the upper surface (like the fluid ve- 
locity) is also expressed in terms of the isothermal sound speed. 
With these scalings, the governing equations become 
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Table 1. The non-dimensional parameters 



Parameter 


Definition 


Values used 


r 


Cp/Cv 


5/3 


m 


gd/R,AT - 1 


1.0 





at/t 


4.0 


K 


K/p dcp(R,T ) [ > 2 


0.0438 


£> 


rjpoCp/K 


0.2, 0.1 or 0.05 


a 


fxcp/K 


0.1 



+ !<2LJ2( (rT 2 /2 + ^|vxB| 2 ), 



<9B 

— = Vx(uxB - < V x B) , V • B = 0. 

at 



(4) 



The pressure, P, satisfies P = pT, whilst the components of the 
stress tensor, r are given by 



an, on : 2 i)ii h ^ 
,J dxj dxj 3 dxk ' J ' 



Note that these equations have a simple equilibrium solution, 
corresponding to a hydrostatic polytropic layer, in which u = 
B = 0, T = 1 + #z and p = T m . 



2.2. Model parameters 

Various non-dimensional parameters appear in this model. These 
are defined in Table [TJ Here, y — 5/3 is the ratio of the specific 
heat capacities, m - 1 .0 is the polytropic index, whilst 6 = 4.0 
is a measure of the thermal stratification. This choice of param- 
eters implies that the layer is superadiabatically-stratified, with 
the temperature varying by a factor of 5 across the layer. The pa- 
rameter a- is a non-dimensionalised thermal diffusivity, <x is the 
Prandtl number, whilst £o represents the ratio of the magnetic to 
the thermal diffusivity at the top of the layer. To interpret these 
coefficients, it is convenient to introduce the Reynolds number 



<Re 



Pmid U n 



Kcr 



(6) 



(where U rms is the rms velocity and p„„j is the horizontally- 
averaged mid-layer density) and the magnetic Reynolds number 



Km = 



(7) 



Note that the depth of the layer, which equals unity in these di- 
mensionless units, has been used as the characteristic lengthscale 
in the definitions of these Reynolds numbers. With the parame- 
ter values that are given in Table [TJ He « 420 for the fully- 
developed hydrodynamic flow, which implies that the convection 
is highly turbulent. We have chosen three values for f > which 
imply that <Rm « 70 (for ft = 0.2), <Rm ~ 140 (for ft = 0.1) or 
Rm « 280 (for ft = 0.05). Note that the magnetic Prandtl num- 
ber, Pm = fbnl'Re is either 0.17, 0.33 or 0.67. As described in 
the Introduction, Pm <K 1 throughout the solar convection zone. 
Although that parameter regime cannot be reached in direct nu- 
merical simulations, we have at least ensured that Pm < 1 in 
all cases. The exploration of the "low" Pm regime is one of the 
novel aspects of this study. 
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Fig. 1. The temperature distribution in a horizontal plane just be- 
low the upper surface of the domain, at t — 0. Brighter contours 
(5) correspond to regions of warmer fluid. 



2.3. Initial conditions 

All of the calculations that are considered in this paper are based 
upon a particular hydrodynamic flow (as determined by the pa- 
rameters that are given in Table [T]). This is evolved in time until 
the convection has reached a statistically-steady state. A snap- 
shot of the resulting flow is shown in Figure [TJ This plot shows 
the temperature distribution in a horizontal plane just below the 
upper surface of the computational domain. Bright contours cor- 
respond to warm upflows, whilst dark contours correspond to 
cooler convective downflows. The flow is characterised by a 
time-dependent pattern of granular convection, reminiscent of 
the one observed at the solar surface. Once the convection has 
reached a statistically-steady state, we insert a magnetic field 
into the flow. In what follows, we define this insertion time to 
be t = 0. This magnetic field takes the form of a weakly twisted, 
horizontal magnetic flux tube that is aligned with the y-axis. The 
radius of the tube, R, is assumed to be 0.1 (which is 10% of the 
depth of the domain) in most cases, although a small number of 
calculations are carried out with a wider tube (R = 0.15) . The 
axis of the flux tube is centred at x = 4.0 and z = 0.75. The three 
components of this field are given by 



B y = BoeH'V*) 2 ] 
B x = 0.4(z - 0.75)B V 
B z = -0.4(x - 4.0)By. 



(8) 



where r = y/(x - 4.0) 2 + (z - 0.75 ) 2 and B is a free parameter 
that can be adjusted to control the strength of the magnetic field. 
The factor of 0.4 in the definitions of B x and B z implies that the 
initial flux tube is weakly twisted. The effect of increasing the 
twist of the tube (by increasing this factor) is briefly discussed 
in the next section. 

The strength of the magnetic field is a free parameter at this 
stage. Recalling that we are working with dimensionless vari- 
ables, it is important to relate the value of Bo to the properties of 
the flow. We define the plasma B to be the ratio of the gas pres- 
sure in the unperturbed polytrope at the initial depth of the tube 
(i.e. at z = 0.75) to the magnetic pressure. The three cases that 
are considered in detail in this paper are Bq = 0.96, Bq = 2.53 
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and Bo = 4.62, which correspond to B — 35, B — 5 and B — 1.5 
respectively. It is also of interest to relate the initial magnetic 
energy density of the field to the kinetic energy density of the 
convection. At a depth of z = 0.75, a value of Bq = 1.25 would 
give a peak field strength that would be in approximate equipar- 
tition with the kinetic energy of the local convective motions. 
However, given that the motion is dominated by strong convec- 
tive downflows, it makes more sense to define an equipartition 
field strength based upon the mean kinetic energy of the convec- 
tive downflows at this depth. This implies an equipartition value 
of B eq = 0.96, which (in turn) implies that the three cases that 
are considered in detail in this paper correspond to Bo ~ B eq , 
Bq sa 2.6B eq and Bo « 4.8B e? . In this sense, this range of values 
for Z?o is similar to those considered in previous studies. 

Magnetic fields of this strength will exert a dynamically- 
significant Lorentz force upon the flow. For a weakly twisted 
flux tube of this form, the x and z components of the magnetic 
field are much smaller than the B y component. As a result of 
this, the dominant component of the Lorentz force is due to ra- 
dial gradients in the magnetic pressure which tend to force fluid 
away from the tube axis. To compensate for this magnetic pres- 
sure, we reduce the gas pressure inside the magnetic flux tube 
at t — in such a way that the total (gas plus magnetic) pres- 
sure distribution inside the tube matches the gas pressure distri- 
bution just before the field was introduced. Although the tube 
is not in perfect force balance with its surroundings, this sim- 
ple metho d produces a very g ood approximation to such a state. 
Following Cheung et al. (2007) we also modify the local veloc- 
ity field within the flux tube so as to minimise convective per- 
turbations to the flux tube during the early stages of evolution. 
Using the notation that we adopted for the magnetic field, we set 

«„ = (l-eB* 2 ])„. (9) 

As described in ICheung et al.l d2007h . this corresponds to a ve- 
locity field that vanishes along the axis of the tube, smoothly 
matching up to the original velocity field near the edge of the 
tube. 

Having specified the gas pressure inside the flux tube, we are 
free to specify one other thermodynamic variable. When consid- 
ering buoyancy instabilities, it is natural to discuss the entropy 
of the flux tube. We specify an initial entropy distribution S„ ew 
that is maximal (and uniform) along the tube's axis, but matches 
smoothly onto the background entropy profile, S , at its outer ra- 
dius. Hence we set 

S„ ew = s(l-e[- ( ^]) + 5oeH^]. (10) 

Normalising the entropy by S mean (0.75), the horizontally- 
averaged entropy at z = 0.75, we consider two different cases. 
In the lower entropy case, the peak entropy along the axis of 
the tube is given by So = 1.26S mean (0.15). This value of Sq is 
comparable to the peak entropy in the convective upflows at this 
depth. In the higher entropy case, we choose a value for So that 
is 50% higher than this, i.e. So — 1.9S mea „(0.75). Note that we 
have also considered some other cases (not shown in this pa- 
per) in which the entropy was uniform in the region r < 0.1. 
In fact, these simulations evolved in a very similar way to the 
radially varying cases. This is presumably due to the fact that 
thermal diffusion rapidly smooths out steep gradients in temper- 
ature (which, in turn, smooths out any sharp radial gradients in 
the entropy distribution). Given that the evolution of these simu- 
lations seems to be rather insensitive to the precise radial distri- 
bution of entropy in the vicinity of the flux tube, only the radially 
non-uniform cases will be discussed in this paper. 



Table 2. A summary of the simulations 



Simulation 


S q/S mean 


R 


P 


Re 


<Rm 


Pm 


LI 


1.26 


0.1 


1.5 


420 


140 


0.33 


L2 


1.26 


0.1 


5.0 


420 


140 


0.33 


L3 


1.26 


0.1 


35.0 


420 


140 


0.33 


L4 


1.26 


0.1 


1.5 


420 


70 


0.17 


L5 


1.26 


0.1 


5.0 


420 


70 


0.17 


L6 


1.26 


0.1 


35.0 


420 


70 


0.17 


L7 


1.26 


0.1 


1.5 


420 


280 


0.67 


L8 


1.26 


0.1 


5.0 


420 


280 


0.67 


L9 


1.26 


0.1 


35.0 


420 


280 


0.67 


WT1 


1.26 


0.15 


1.5 


420 


140 


0.33 


WT2 


1.26 


0.15 


1.5 


420 


70 


0.17 


WT3 


1.26 


0.15 


1.5 


420 


280 


0.67 


HI 


1.9 


0.1 


1.5 


420 


140 


0.33 


H2 


1.9 


0.1 


5.0 


420 


140 


0.33 


H3 


1.9 


0.1 


35.0 


420 


140 


0.33 


El 




0.1 


1.5 


420 


140 


0.33 


E2 




0.15 


1.5 


420 


140 


0.33 



Table [2] summarises the numerical simulations in this study. 
Most of the simulations (L1-L9) correspond to the lower en- 
tropy, thin tube (R = 0.1) cases, for different values of B and 
different values of Rm. HI, H2 and H3 illustrate the effects of 
increasing the initial peak entropy within the thin flux tube (but 
are otherwise identical to cases LI, L2 and L3). WT1, WT2 
and WT3 represent three calculations that were carried out with 
a wider flux tube, using the lower entropy initial conditions. 
Finally, in the El and E2 simulations, the velocity field, the gas 
pressure and the entropy distribution are exactly as they were 
before the field was introduced. Thus, no modification is made 
to the convection at t = 0. This implies that the Lorentz force 
will significantly perturb the flow during the early stages of evo- 
lution. The rationale for these particular calculations will be dis- 
cussed later in the paper. 

We illustrate the initial condition for one of the low entropy 
cases (LI) in Figure [2] which shows the distribution of various 
quantities in the x = 4.0 plane. The upper plot shows contours 
of constant B y . The localised nature of this initial flux tube is 
clearly apparent. The middle plot in Figure [2] shows contours of 
constant entropy. It is clear from this plot that the peak entropy in 
the flux tube is comparable to the peak entropy in the convective 
upflows. We would expect a tube with this entropy distribution 
to be buoyant relative to its surroundings. Finally, the lower part 
of Figure [2] shows the density perturbation in this plane. To ob- 
tain this plot, we have subtracted off the polytropic background 
stratification in order to highlight the variations in p due to the 
magnetic flux tube. As expected, the density perturbation in this 
strong field (B — 1.5) case takes its minimum value within the 
flux tube. Here, the minimum value of p is approximately 40% 
of the mean density at this depth. The associated reduction in 
density is smaller in the B = 5 case, and almost non-existent in 
the B = 35 case. This simply reflects the fact that the magnetic 
field is weaker in the higher B cases. 

To evolve these simulations, we use a well tested mixed 
pseudospectral/finite difference code. All horizontal derivatives 
are evaluated in Fourier space, using standard Fast Fourier 
Transform (FFT) libraries. Fourth-order finite differences are 
used to calculate vertical derivatives. The time-stepping is car- 
ried out using a third-order Adams-Bashforth scheme, with a 
variable time-step. The code is parallelised using MPI. All sim- 
ulations are carried out ona512x512x96 mesh. 
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Fig. 2. Vertical slices through the x = 4.0 plane, showing con- 
tours of B v (top), entropy (middle) and the density perturbation 
(bottom) for the LI case at t — 0. 



3. Results 

3.1. Varying the initial conditions 

We first consider cases LI, L2 and L3, which correspond to the 
lower entropy simulations at 'Rm ss 140, for B — 1.5, B = 5 and 
B — 35 respectively. In all of these cases, convective motions 
play a significant role in determining the early evolution of the 
system, disrupting the flux tube by pulling out loops of magnetic 
field that are then advected around the domain. The early stages 
of convective disruption are illustrated in Figure [3] which shows 
contours of B y /Bo in the x = 4.0 plane for each of these cases 
at t — 0.36. Normalising the magnetic field by Bq allows us to 
compare cases with different values of B. Defining the convec- 
tive turnover time, T conv , to be the inverse of the rms velocity, 
T com , ~ 1.56 in these dimensionless units. Hence t = 0.36 cor- 
responds to a small fraction of the turnover time. However, even 
after this short time period, we see that the flux tube has been 
perturbed in a significant way by the convective flows. This con- 
vective disruption is most apparent in the weaker field (higher B) 
cases, where small-scale perturbations to the magnetic field have 
produced a magnetic field with an undulating sea-s erpent config- 
uratio n (as observed in previous calculations, e.g. lCheung et al.l 
120071) . Only in the B — 1 .5 case does the field appear to be strong 
enough to partially resist the convective disruption on this time- 
scale. Obviously stronger fields possess a higher level of mag- 
netic tension, which makes it more difficult for the surrounding 
convective motions to perturb the magnetic flux tube. 

The subsequent evolution of cases LI and L3 is illustrated 
in Figure 3] which shows the contours of B v /Bq in the x = 4.0 
plane, for each case, at t — 0.80 and t = 1.59. At these later 
times, the differences between the stronger field and the weaker 
field case have become more pronounced. Focusing initially 
upon the weaker field case (L3), it is clear that the flux tube 
has been further disrupted by the convective motions. In the 
t = 0.80 plot, the upflows and downflows in the lower part of 
the domain have produced many loop-like structures, with sev- 
eral well-defined regions of negative B y (which correspond to the 
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Fig. 3. Vertical slices through the x = 4.0 plane, showing con- 
tours of By/Bo at t — 0.36 for three different cases, LI (B = 1.5, 
top), L2 (B = 5.0, middle) and L3 (B — 35, bottom). In each case, 



the contours are evenly spaced in the range 
0.9. 



-0.13 < B Y / B < 



dark contours in Figure 0]). Given the relatively weak magnetic 
fields that are present at this stage of the simulation (even the 
peak horizontal field strength is now well below the equipartition 
value that was discussed in the previous section), we would ex- 
pect the magnetic field to play a relatively passive role in the dy- 
namics. The efficiency with which the convective motions have 
advected the magnetic flux around the domain would appear to 
be consistent with this idea. We can see from the t = 1.59 plot 
for the L3 case that a strong upflow in the vicinity of y = 5 has 
advected some horizontal magnetic flux into the upper layers of 
the computational domain. Although convective disruption does 
also play a role in the stronger field case (LI), there are fewer 
loop-like structures in the magnetic field distribution at t = 0.80, 
which suggests that the magnetic flux tube maintains a greater 
degree of coherence than in the weaker field case. Furthermore, 
as illustrated in the t = 1.59 plot, a much broader concentra- 
tion of magnetic flux rises towards the surface regions in this 
case than it does in the corresponding weaker field calculation. 
Even though the peak field strength in this emerging flux con- 
centration is a relatively small fraction of its initial value, Bo, it 
is still comparable to the equipartition value, which suggests that 
it should still exert a dynamically-significant Lorentz force upon 
the flow. We shall discuss the flux emergence pattern for the LI 
case in more detail in the next subsection. 

Figures [3] and [4] clearly show that there is some dependence 
upon B in these simulations. The evolution of the flux tube 
in the weakest field case (L3) appears to be largely (possibly 
completely) determined by the surrounding convective motions. 
However, further quantitative analysis is needed in order to de- 
termine whether or not magnetic buoyancy is influencing the 
evolution of the flux tube in the lower B cases. In order to investi- 
gate this possibility, we have traced the evolution of the unsigned 
emerging magnetic flux at z — 0. This is given by 



O(f) 



Jo Jo 



\B z (z = Q,t)\dxdy. 



(11) 
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Fig. 4. Contours of B v /B in the x - 4.0 plane. Top: LI (J3 = 1.5) at t = 0.80 (left) and t - 1.59 (right). Bottom: The same plots for 
the L3 (B = 35) case. In the t = 0.80 plots, the contours are evenly spaced in the range -0.25 < B y /Bo < 0.45, whilst the contour 
range is -0.22 < B y /B < 0.26 in the t = 1.59 plots. 
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Fig. 5. The magnitude of the total unsigned magnetic flux, <J>(f), 
at z = 0, plotted against time for simulations LI, L2 and L3 
(the lower entropy, thin tube, *Rm ss 140 cases). The solid line 
corresponds to LI (B = 1.5), the dotted line corresponds to L2 
(B = 5), whilst the dashed line corresponds to L3 (B — 35). In 
each case, the flux has been normalised by the initial magnetic 
flux in the y-direction. 



So far only the strength of the imposed field has been varied. 
As indicated in Table|2] cases HI, H2 and H3 adopt the same val- 
ues of B and fbn as LI, L2 and L3 (respectively), but this time 
with a higher peak entropy within the initial flux tube. Figure [6] 
shows the flux emergence as a function of time for these cases. 
Comparing this plot with Figure|5] it is clear that there is gener- 
ally an enhanced rate of flux emergence in these higher entropy 
calculations. This is to be expected given that the higher entropy 
flux tube should be more buoyant than those considered in cases 
LI, L2 and L3. However, it is again notable that there is only 
a modest dependence upon B in the weaker field cases, H2 and 
H3, which exhibit very similar flux emergence rates. As in the 
lower entropy case, there is again a slightly higher rate of flux 
emergence in the HI (J3 = 1.5) case. As before, this suggests 
that magnetic buoyancy plays a role in the stronger field calcula- 
tion. We can conclude from this that the enhancement of the flux 
emergence rate in the stronger field case is a robust result that 
does not depend critically upon the choice of the initial entropy 
distribution along the axis of the flux tube. Some simulations 
with even higher values of the initial entropy were considered. 
However the emerging flux in these cases tended to be associ- 
ated with regions of exceptionally high temperature (something 
that is not observed in the Sun), so these calculations were not 
pursued in any detail. 



3.2. Varying the flux tube geometry 



To compare cases at different values of B, we normalise this 
quantity by the initial flux in the y-direction. The results of this 
quantitative analysis for LI, L2 and L3 are shown in Figure [5] 
Cases L2 and L3 both give very similar flux emergence rates, 
which indicates that the flux emergence in each case is driven 
almost entirely by the convective upflows. The marginal differ- 
ence between the flux emergence rates in these two cases would 
then simply reflect the fact that stronger fields are more resistant 
to advection by the convective motions. However, the slightly 
higher flux emergence rates in simulation LI can only be ex- 
plained if magnetic buoyancy is enhancing the flux emergence 
rate in this case. We therefore conclude that magnetic buoyancy 
is playing a role in the evolution of the flux tube in this lower 
B (stronger field) calculation, helping the convective upflows to 
advect dynamically-significant concentrations of magnetic flux 
into the upper layers of the domain. 



We also considered variations in the twist and the radius of 
the initial magnetic flux tube. Our aim here was to determine 
whether it was possible to change the geometry of the tube in 
such a way as to enhance the effects of magnetic buoyancy. 
Given that this was the only lower entropy case in which there 
are any indications of magnetic buoyancy, we used the strong 
field simulation LI (B — 1.5) as the basis for investigating these 
variations in the flux tube geometry. In this case, increasing the 
original twist by a factor of 10 does not seem to influence the 
evolution of the tube, or the flux emergence rate. Although this 
enhancement of the twist does imply that the fieldlines now 
make approximately 3-4 full turns around the axis of the im- 
posed flux tube, the resulting magnetic tension is still much 
smaller than the initial magnetic pressure gradient across the ra- 
dius of the tube. It is therefore not surprising that the dynamics 
are largely unaffected by this change. 
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Fig. 6. As figure [5] but this time for the high entropy cases at Fig. 8. As figure [5] but for simulations WT1 and LI. The solid 

Rm » 140 (HI, H2 and H3). The solid line corresponds to HI line corresponds to the wide tube case (WT1), whilst the dotted 

(B = 1.5), the dotted line corresponds to H2 (B = 5), whilst the line corresponds to the thin tube case (LI), 
dashed line corresponds to H3 (B — 35). 
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Fig. 7. Vertical slices through the x = 4.0 plane for the WT1 
simulation, showing the density perturbation (top) and contours 
of constant B v /Bq (middle) at t — 0.36. The lower plot shows 
contours of constant B v /Bo at t = 0.80. In the middle plot, the 
contours are evenly spaced in the range -0.09 < B y /Bo < 0.85. 
In the lower plot, the corresponding range is -0.06 < B v /Bq < 
0.55 

Increasing the width of the flux tube has a much more sig- 
nificant effect upon this system. Case WT1 is identical to LI in 
every respect except for the value of R, which has been increased 
from 0.1 to 0.15. This more than doubles the initial horizon- 
tal magnetic flux in the initial condition. The upper two plots in 
Figure|7]show the density perturbation and the contours of By/Bg 
in the x = 4.0 plane at t — 0.36. At this instant in time there is (as 
in case LI) some indication of modest convective disruption, al- 
though there is still a coherent flux tube with a pronounced den- 



sity deficit. Indeed, at z — 0.75, the minimum density is approx- 
imately 60% of the mean density at this depth. So even after the 
initial phase of evolution, this is still a configuration that should 
be susceptible to magnetic buoyancy. The lower plot in Figure [7] 
shows the distribution of B v /Bq in the same plane at t — 0.8. This 
plot shows that some coherent flux concentrations have risen into 
the upper layers of the domain. Comparing this plot with the cor- 
responding plot for LI at t — 0.80 (which is shown in the top left 
part of Figure |4]i, we see that there seems to be more flux in the 
upper layers of the domain in the wider tube case. To illustrate 
this point in more quantitative terms, Figure [8] shows the (un- 
signed) flux emergence as a function of time for the WT1 simu- 
lation, plotted alongside the equivalent flux emergence curve for 
LI. In each case, the flux is normalised by the initial imposed 
flux in the y-direction. Given that the (normalised) flux emer- 
gence rate is significantly higher in case WT1, we can conclude 
that magnetic buoyancy is playing a much more significant role 
in the evolution of this wider magnetic flux tube. This suggests 
that the efficiency of magnetic buoyancy relative to the effects of 
convective disruption is determined not only by the initial peak 
magnetic field strength, but also by the width of the imposed flux 
tube. 

The left-hand side of Figure|9]shows the flux emergence pat- 
terns for cases LI and WT1 at t = 1.59. In the LI case, the 
vertical magnetic field distribution in the near-surface region is 
dominated by two well-defined bipolar magnetic regions. As il- 
lustrated in Figure [9] these flux concentrations emerge initially 
within the warm granular interiors, in regions of convective up- 
flow. However once the magnetic fields reach the near-surface 
layers, the horizontal convective motions start to sweep these 
flux concentrations into the intergranular lanes. So at later times 
these emerging vertical flux concentrations accumulate in the 
convective downflows. The flux emergence pattern in case WT1 
is similar, although more flux has emerged, and the resulting 
magnetic field distribution is much less localised than the flux 
emergence pattern from case LI. The volume renderings on the 
right-hand side of Figure [9] show isosurfaces of the magnetic 
field strength (|B| = 0.25). Although the flux tube is wider in 
the WT1 case, the flux tube structure is similar in both cases. It 
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Fig. 9. Simulations LI (top) and WT1 (bottom) at t — 1.59. The plots on the left-hand side show the temperature distribution in 
a horizontal plane just below the upper surface of the computational domain (shaded contours) and the distribution of B, in the 
same horizontal plane (black and white contours). The volume renderings on the right-hand side of this figure show isosurfaces of 
|B| = 0.25. Some fieldlines have also been traced out along these isosurfaces in order to aid visualisation. 



is clear that the emerging bipoles that are observed at the surface 
belong to the same undulating flux tube, and are therefore con- 
nected by subsurface magnetic field lines. So in each case, the 
flux tube has resisted convective disruption effectively enough 
to maintain some semblance of its initial structure throughout 
the flux emergence process. 

3.3. Varying the magnetic Reynolds number 

All of the calculations that have been described so far have a 
magnetic Reynolds number of Rm « 140. Although the ohmic 
decay time based on the depth of the layer, r n , is approximately 
140 convective turnover times for Rm m 140, the decay time 
based upon the width of the flux tube is approximately two or- 
ders of magnitude smaller than that. Hence, magnetic dissipation 
must be playing some role in the evolution of the system on the 
timescales of interest. The crucial question is whether or not the 
evolution of the system is dominated by magnetic dissipation. 
To investigate this issue, we varied the value of Rm. Cases L4, 
L5 and L6 are identical to LI, L2 and L3 in every respect except 
for the fact that Rm w 70, as opposed to Rm » 140. In cases L7, 
L8 and L9, we adopted a higher value of Rm « 280. We also 
carried out some wider tube calculations (WT2 and WT3) with 
these values for the magnetic Reynolds number. As described 
in the Section 2.2, these choices for Rm ensure that Pm < 1 in 



all cases, and this restriction upon Pm is the main reason why 
higher values of Rm were not considered in the present study. 

Figures[10]and[TT]show the normalised flux emergence rates 
for the thin flux tube, lower entropy cases with Rm as 70 and 
Rm ~ 280 respectively. These plots should be compared with 
Figure|5] which shows the corresponding plot for the Rm as 140 
cases. The first point to notice is that that B dependence of the 
flux emergence rate is strongly-dependent upon Rm. In partic- 
ular, the enhancement of the flux emergence rate due to mag- 
netic buoyancy in the B = 1.5 simulations becomes much more 
pronounced at higher values of the magnetic Reynolds number. 
After one convective turnover time (at t as 1.55), the normalised 
surface fluxes in the Rm as 70 simulations differ by only a few 
percent whereas in the Rm as 280 simulations the normalised 
fluxes differ by approximately 15% after a similar time period. 
We can conclude from this that magnetic buoyancy is more ef- 
ficient at higher magnetic Reynolds number (presumably a con- 
sequence of the fact that there is less magnetic dissipation), pro- 
vided that the magnetic field is strong enough initially to pro- 
duce a magnetic pressure that is a significant fraction of the lo- 
cal gas pressure. However it is not just the efficiency of magnetic 
buoyancy that is dependent upon the magnetic Reynolds number. 
Comparing Figures [10] and [TTJ we see that the flux emergence 
rates in the higher magnetic Reynolds number cases are system- 
atically lower than the corresponding flux emergence rates in the 
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Fig. 10. As figure[5] but this time for the Rm a 70 thin tube cases 
(L4, L5 and L6). The solid line corresponds to L4 (J3 = 1.5), the 
dotted line corresponds to L5 (J3 = 5), whilst the dashed line 
corresponds to L6 ifi = 35). 
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Fig. 11. As figure |5j but this time for the Rm a 280 thin tube 
cases (L7, L8 and L9). The solid line corresponds to L7 (ft = 
1 .5), the dotted line corresponds to L8 (J3 — 5), whilst the dashed 
line corresponds to L9 (J3 = 35). 



smaller Rm calculations. Although not shown here, the same is 
true in the wider tube calculations. The reason for this reduced 
flux emergence rate is that convective disruption (like magnetic 
buoyancy) is also more efficient at higher magnetic Reynolds 
numbers, because there is a greater tendency for magnetic field 
lines to be advected by the flow. As a result of this, loops of flux 
are continually removed from the tube, before being reprocessed 
by the local convective motions. This causes less flux to emerge 
at the surface of the domain. As described in the Introduction, we 
know from previous studies (lDorchetal.1 12001) that convective 
motions generally reduce the efficiency of flux emergence. We 
can now make the stronger statement that this reduction in the 
flux emergence rate due to convective disruption is more pro- 



nounced at higher magnetic Reynolds numbers in the Pm < 1 
parameter regime. 

Given that the flux emergence rates are strongly dependent 
upon Rm, we would also expect the flux emergence patterns to 
show a similar 7?m-dependence. Figure [12] shows the flux emer- 
gence patterns at t — 1 .59, alongside a volume rendering of the 
associated magnetic flux tubes, for cases L7 and WT3 (the low- 
est p, lower entropy cases for Rm as 280). Although the strongest 
emerging flux concentrations are still found in the vicinity of 
x = 4.0, which corresponds to the location of the axis of the ini- 
tial magnetic flux tube, many other weaker magnetic flux con- 
centrations are also observed. This indicates that the flux tube 
has undergone more convective shredding in these higher mag- 
netic Reynolds number cases. This convective disruption is high- 
lighted even more clearly in the volume renderings, particularly 
in the wide tube case where the magnetic field distribution is 
highly disordered and fragmented. Although not shown here, 
very little convective shredding is observed in the Rm « 70 
cases, and the flux emergence pattern is dominated by a small 
number of coherent dipoles. 

3.4. Non-equilibrium expanding flux tubes 

In this subsection, we consider the last two cases listed in 
Table |2j namely El and E2. Here, we make no initial pertur- 
bation to the flow, pressure or entropy at t — 0. The idea of these 
simulations is to allow the convection to select its own "initial 
condition", by adjusting itself in response to the presence of a 
dynamically-significant magnetic field. Any magnetic buoyancy 
is therefore induced in a self-consistent way, albeit from a non- 
equilibrium initial condition. This approach, whilst clearly over- 
simplified (and unrealistic for the solar convection zone), has 
the great advantage that it is not necessary to specify a partic- 
ular distribution for any of the thermodynamic quantities along 
the tube, in contrast to the other cases that are considered in this 
paper. Given that there is considerable uncertainty regarding the 
most suitable initial conditions for idealised calculations of this 
type, we would argue that simulations El and E2 are useful il- 
lustrative calculations in this context. 

Obviously, the initial phase of evolution of these simulations 
is dominated by the effects of the Lorentz force. The radial mag- 
netic pressure gradient tends to drive fluid away from the axis 
of the tube. Furthermore, this occurs on an Alfvenic timescale, 
which means that this process occurs very rapidly (in these cases, 
the Alfven speed at the depth of the flux tube is approximately 4 
times larger than a typical convective velocity). This pressure- 
driven radial outflow partially evacuates the flux tube. It also 
causes the flux tube to expand in the radial direction. By con- 
servation of flux, this expansion leads to a reduction in the mean 
magnetic field strength within the magnetic flux tube. As a result 
of this expansion, the flux tube rapidly reaches a state in which 
it is in (near) total pressure balance with its surroundings. Given 
that this adjustment is very rapid, convective perturbations play 
a negligible role during this early phase of evolution. 

Having established that these cases rapidly evolve towards 
pressure-balanced states, we now consider the subsequent evo- 
lution of each of these calculations. The initial magnetic field in 
case El is identical to that of case LI, with fi — 1.5 and R =0.1. 
The top two plots in Figure Q~3] show a snapshot of this case at 
t = 0.36. The upper plot shows the density perturbation in the 
x = 4.0 plane, whilst the middle plot shows the B y distribution 
in the same plane (normalised by Bq). At this instant in time, 
the minimum density within the flux tube is approximately 75% 
of the mean density at z = 0.75. Although this level of partial 
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Fig. 12. As figure[9] but here for the corresponding Rm a 280, lower entropy cases, L7 (top) and WT3 (bottom). 



evacuation is not negligible, it is modest compared to some of 
the other cases. The B y /Bo distribution at t = 0.36 is compara- 
ble to that shown in the upper plot of Figure [3] (where the same 
contour spacings were used). The lower plot in Figure [T3lshows 
the By/Bo distribution in the x = 4.0 plane at t = 0.8. Again this 
horizontal magnetic field distribution is comparable to that of the 
corresponding plot for the LI case (shown in the top left plot in 
Figure |4]i. In qualitative terms, these plots clearly indicate that 
the flux tube in this new calculation has evolved in a rather sim- 
ilar way to the flux tube from the LI case, despite the unrealistic 
initial conditions. Figure [14] shows the corresponding results for 
the E2 case (B — 1.5 and R = 0.15), which can be compared 
directly to WTl (as illustrated in Figure|7]i. At t — 0.36, the min- 
imum density at z — 0.75 for this E2 case is approximately 70% 
of the mean density at this depth, which is slightly higher than 
the minimum density at this level at the corresponding stage of 
the WTl simulation. However, in almost every other respect, the 
results from this E2 case are qualitatively rather similar to those 
of the WTl simulation. The flux emergence rates for cases El 
and E2 are shown in Figure[l5] These are clearly comparable to 
the corresponding flux emergence curves for the LI and WTl 
cases. So although the initial conditions in cases El and E2 are 
rather idealised (and probably unrealistic), we have shown that 
these simulations do evolve in a similar way to our previous cal- 
culations. This suggests that local simulations of this type are 
relatively insensitive to the precise choice of initial conditions, 
with the evolution of the system (for a given value of Rm) being 
determined primarily by the initial magnetic field distribution. 



Given the wide range of plausible choices of initial conditions in 
local models of this type, this is an encouraging result. 



4. Summary and discussion 

We have used numerical simulations to investigate the dynami- 
cal evolution of a magnetic flux tube embedded within a granular 
convective layer. Throughout the calculations we have adopted 
the same basic hydrodynamic flow, with a Reynolds number of 
Re » 420. When the magnetic flux tube is introduced, the gas 
pressure is adjusted so that the tube is in pressure balance with 
its surroundings. The tube is then made buoyant by choosing 
an appropriate entropy distribution. The subsequent evolution of 
this tube is influenced by two key processes, namely convective 
disruption and magnetic buoyancy. However the extent to which 
magnetic buoyancy plays a role in the flux emergence process 
depends in a rather subtle way upon the model parameters. In 
this systematic survey, we varied the initial peak magnetic field 
strength, the twist and width of the flux tube, the magnitude of 
the entropy along the tube axis and the magnetic Reynolds num- 
ber. In all cases, the magnetic Prandtl number, Pm - RmIRe 
is less than unity, as would be expected for the solar convection 
zone. This is the first study of this type to explore this low Pm 
regime. 

At a moderate magnetic Reynolds number of Rm » 140, the 
evolution of a thin magnetic flux tube is partially determined by 
the effects of convective disruption. In the weak field (high B) 
cases, the magnetic field appears to play a relatively passive role 
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Fig. 13. Vertical slices through the x = 4.0 plane for the El sim- 
ulation, showing the density perturbation (top) and contours of 
constant B y /Bo (middle) at t — 0.36. The lower plot shows con- 
tours of constant B y /Bo at t — 0.80. In the middle plot, the con- 
tours are evenly spaced in the range -0.13 < B y /Bo < 0.9. In the 
lower plot, the corresponding range is -0.25 < B y /Bo < 0.45. 





Fig. 14. Like Figure [13] but this time for the E2 simulation. In 
the middle plot, the contours are evenly spaced in the range 
-0.09 < By/Bo < 0.85. In the lower plot, the corresponding 
range is -0.06 < B y /B < 0.55. 



in the dynamics. However, there is an enhanced flux emergence 
rate in the strongest field cases due to the effects of magnetic 
buoyancy. Increasing the magnitude of the initial peak entropy 
along the tube axis does not change the B dependence, but it 
does increase the rate of flux emergence in all cases. Varying the 
twist of the flux tube does not seem to influence the evolution 
of the system (presumably because the magnetic tension is still 
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Fig. 15. As figure [5] but for simulations El and E2. The solid 
line corresponds to the wide tube case (E2), whilst the dotted 
line corresponds to the thin tube case (El). 



much weaker than the magnetic pressure gradients in the flux 
tube). However, increasing the width of the tube does produce 
an enhanced flux emergence rate in the strong field case. This is 
presumably because a thicker tube is more robust, and tends to 
resist convective disruption more effectively than a thinner flux 
tube. As a result of this, the flux tube maintains its coherence 
more effectively during the early stages of evolution, which al- 
lows magnetic buoyancy to operate in a more efficient manner. 
So our main conclusion from the %n m 140 cases is that the con- 
tribution of magnetic buoyancy to the flux emergence process is 
directly related not only to the peak field strength, but also to the 
width of the initial magnetic flux tube. 

We also varied the magnetic Reynolds number in this system 
in order to determine the extent to which the evolution of the tube 
is influenced by the effects of magnetic dissipation. This para- 
metric survey produced two key results. Firstly, magnetic buoy- 
ancy contributes most effectively to the flux emergence process 
in the low B, high magnetic Reynolds number regime. Secondly, 
for a given value of B, the flux emergence rate is always lower in 
the higher Km regime. This is due to the fact that convective dis- 
ruption is much more efficient at higher magnetic Reynolds num- 
ber (where field lines are more easily advected by the flow). The 
extent to which this result depends upon the magnetic Prandtl 
number is unclear, although moving to higher Km (thereby in- 
creasing Pm) would be computationally simpler than reducing 
Pm by increasing "Re. At higher fbn this convective flow should 
be capable of sustaining a small-scale dynamo, although it is not 
clear how such a flow would interact with an emerging flux tube. 
This is a possible area for future work. 

Motivated by some of the uncertainty surrounding the most 
appropriate choice of initial conditions for problems of this type, 
we also carried out some idealised simulations in which the gas 
pressure was not adjusted when the flux tube was introduced. 
We stress that this is not a realistic situation for the solar con- 
vection zone given that this leads to a strong imbalance between 
the magnetic pressure and the local gas pressure. However, this 
does allow the system to adjust to the presence of the flux tube 
in a self-consistent way. This adjustment occurs very rapidly and 
(somewhat remarkably) the subsequent evolution of the system 
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is comparable to that of the previous set of simulations, at least 
at comparable magnetic Reynolds numbers. Results from these 
idealised calculations strongly suggest that the only aspect of the 
initial conditions that influences the evolution of the flux tube 
is the initial magnetic field distribution. This is an encourag- 
ing result because it suggests that flux emergence simulations 
are largely insensitive to the precise choice of initial conditions, 
which is a matter of considerable uncertainty in idealised models 
of this type. 
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